Assignment-RECAST Modern Marketing Mix Modeling

Purpose: Revenue vs AdSpend effect

Author: Gerardo Gandara

Contact: gerardo.gandara@gmail.com

Client: andy@getrecast.com

Code created: 2023-10-11

Last updated: 2023-10-11

Source: https://github.com/ggandara13/myrepository

Comment: https://getrecast.com/modern-media-mix-modeling/.

Libraries

These are the libraries to load

pkg <- c('pivottabler', 'tidyverse','tibbletime','anomalize','timetk')
#install.packages(pkg)

library(pivottabler)
library(lubridate)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(tibbletime)
library(anomalize)
library(timetk)
library(forecast)

Read the revenue and AdSpend

We check the date type of the ‘date’ column and change the type to DATE

The CSV data is avaialble in the github repository

setwd("/Users/gerardogandara/Documents/recast/")

## Read file as CSV
df_spend <- read.csv('https://raw.githubusercontent.com/ggandara13/myrepository/master/acme_spend.csv')
df_revenue <- read.csv('https://raw.githubusercontent.com/ggandara13/myrepository/master/acme_revenue.csv')

# Convert char format to DATE format
df_revenue$date <- as.Date(df_revenue$date, format =  "%m/%d/%y")
df_spend$date <- as.Date(df_spend$date, format =  "%m/%d/%y")

# Check the format and content of each dataframe
str(df_revenue)

'data.frame':   1035 obs. of  5 variables:
 $ X              : int  1 2 3 4 5 6 7 8 9 10 ...
 $ date           : Date, format: "2020-01-01" "2020-01-02" "2020-01-03" ...
 $ revenue_dtc    : num  16000 19306 22052 23177 26169 ...
 $ revenue_amazon : num  50000 60332 68914 72428 81778 ...
 $ revenue_walmart: num  34000 41026 46861 49251 55609 ...

str(df_spend)

'data.frame':   10350 obs. of  4 variables:
 $ X      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ date   : Date, format: "2020-01-01" "2020-01-01" "2020-01-01" ...
 $ channel: chr  "facebook_prospecting" "facebook_retargeting" "google_branded_search" "google_nonbranded_search" ...
 $ spend  : num  4154 2679 474 3817 7421 ...

Check there are not NULL values

paste("Total of NULL in Spend:", sum(is.na(df_spend)))

[1] "Total of NULL in Spend: 0"

paste("\nTotal of NULL in Revenue:", sum(is.na(df_revenue)))

[1] "\nTotal of NULL in Revenue: 0"

paste("\n\nSummary\n")

[1] "\n\nSummary\n"

# you can check with summary as well and see the columns
summary(df_spend)

       X              date              channel              spend            Year     
 Min.   :    1   Min.   :2020-01-01   Length:10350       Min.   :    0   Min.   :2020  
 1st Qu.: 2588   1st Qu.:2020-09-15   Class :character   1st Qu.: 1278   1st Qu.:2020  
 Median : 5176   Median :2021-06-01   Mode  :character   Median : 5268   Median :2021  
 Mean   : 5176   Mean   :2021-06-01                      Mean   : 7740   Mean   :2021  
 3rd Qu.: 7763   3rd Qu.:2022-02-15                      3rd Qu.:11649   3rd Qu.:2022  
 Max.   :10350   Max.   :2022-10-31                      Max.   :36401   Max.   :2022  

summary(df_revenue)

       X               date             revenue_dtc     revenue_amazon   revenue_walmart 
 Min.   :   1.0   Min.   :2020-01-01   Min.   :     0   Min.   :     0   Min.   :     0  
 1st Qu.: 259.5   1st Qu.:2020-09-15   1st Qu.: 41643   1st Qu.:130216   1st Qu.: 88547  
 Median : 518.0   Median :2021-06-01   Median : 72051   Median :225161   Median :153109  
 Mean   : 518.0   Mean   :2021-06-01   Mean   : 72741   Mean   :227062   Mean   :154488  
 3rd Qu.: 776.5   3rd Qu.:2022-02-14   3rd Qu.: 97340   3rd Qu.:303352   3rd Qu.:206279  
 Max.   :1035.0   Max.   :2022-10-31   Max.   :193005   Max.   :415952   Max.   :282848  

sum(is.na(df_spend))/nrow(df_spend)*100

[1] 0

sum(is.na(df_revenue))/nrow(df_revenue)*100

[1] 0

Since we have 1,035 rows, we can remove because that represents less than 1% or just replace with zero

Now, let’s replace with zero since it is a time series

df_revenue[is.na(df_revenue)] <- 0
df_spend[is.na(df_spend)] <- 0


cat("Total of NULL in Spend:", sum(is.na(df_spend)))

Total of NULL in Spend: 0

cat("\nTotal of NULL in Revenue:", sum(is.na(df_revenue)))

Total of NULL in Revenue: 0

1- Which channel had the largest increase in spend in 2022 compared to the same date range in 2021?

Using lubridate to filter the year of the date, also use the pivottabler library to group/aggregate it

Answer: online_video

# print the first value of the dataframe
paste("The channel with highest LIFT is:", row.names(df_summary_final)[1] )

[1] "The channel with highest LIFT is: online_video"

2- In terms of total revenue, are there any anomalous days?

Group by day to see the trend and plot

df_revenue_total <- df_revenue %>%
  mutate(total_revenue = ( df_revenue$revenue_dtc + df_revenue$revenue_amazon + df_revenue$revenue_walmart  ))

# Group by sum using R Base aggregate()
agg_df <- aggregate(df_revenue_total$total_revenue, by=list(df_revenue_total$date), FUN=sum)
colnames(agg_df) <- c('date','revenue') 


#Plot sales over 36 months
ggplot(agg_df, aes(x=date, y = revenue)) + geom_line()

paste("Visually we can see some outliers")

[1] "Visually we can see some outliers"

Now, we can use the anomalties library to identify the values

Use time_decompose() to decompose a time series prior to performing anomaly detection with anomalize(). Typically, anomalize() is performed on the “remainder” of the time series decomposition.

The return has three columns: “remainder_l1” (lower limit for anomalies), “remainder_l2” (upper limit for anomalies), and “anomaly” (Yes/No).

We can print the detail of the anomalies, for example the last 4

tail(df_anomalized %>%
  filter(anomaly == 'Yes'),4)

NA

You can also use the timetk package for dynamic plot

agg_df %>% timetk::plot_anomaly_diagnostics(date,revenue, .facet_ncol = 2)

frequency = 7 observations per 1 week
trend = 92 observations per 3 months

3- In which month of the year does Acme tend to make the most revenue?

We can compute month column before we aggregate and sort to get the popular month

Most profitale month in a year is the month number: 3

# Calculate the month column
monthfunction<-function(dataframe, datecolumn) {
    month(dataframe[,datecolumn])
}
agg_df$Month <- monthfunction(agg_df, "date")

# group by month

by_month <- aggregate(agg_df$revenue, by = list(agg_df$Month), FUN = sum)
colnames(by_month) <- c('month','revenue') 



by_month <- by_month[order(by_month$revenue, decreasing = TRUE), ]

# print the first value of the dataframe
print(by_month)

paste("Most profitale month in a year is the month number:", row.names(by_month)[1] )

[1] "Most profitale month in a year is the month number: 3"

4- Does Acme’s marketing spend tend to follow a similar pattern to revenue?

We need to merge both dataframes to have the TS of Revenue using AdSped as regressor Convert the dataframe to TS then plot, then plot the results

Visually we can confirm they have a similar pattern.

#select the date and total_revenue of Acme
df_only_revenue <- select(df_revenue_total, date, total_revenue)

#aggregate spend of all media
df_spend_total <- aggregate(df_spend$spend, by=list(df_spend$date), FUN=sum)
colnames(df_spend_total) <- c('date','spend') 

df_rev_spend <- merge(df_only_revenue, df_spend_total, by = "date")


#Convert dataframe to time series object using the ts() function
acme_ts <- ts(data = df_rev_spend[,c(2,3)])


# Time plot of both variables
autoplot(acme_ts, facets = TRUE)

5- Based on what you see here, could you make the case that marketing spend is causally related to revenue? Why or why not?

Visually we can see the effect has high correlation, and when we calculate the R2 is 0.91 with p=value basically zero. So, it is statiscally correlated.

Is the relation linear? Yes, form the plot above, the relationship is linear.

now we jsut need to verify that the data from each of the 2 variables (revenue, spend) follow a normal distribution.

df_rev_spend

# Shapiro-Wilk normality test for revenue
shapiro.test(df_rev_spend$total_revenue) # => p-value < 2.2e-16

    Shapiro-Wilk normality test

data:  df_rev_spend$total_revenue
W = 0.93373, p-value < 2.2e-16

# Shapiro-Wilk normality test for spend
shapiro.test(df_rev_spend$spend) # => p-value < 2.2e-16

    Shapiro-Wilk normality test

data:  df_rev_spend$spend
W = 0.95587, p-value < 2.2e-16

# revenue 
ggqqplot(df_rev_spend$total_revenue, ylab = "revenue")

# spend
ggqqplot(df_rev_spend$spend, ylab = "spend")

NA
NA

Since the p-value is practically zero, we can reject the null hypothesis that it is normal distributed both variables

so we can use another test to check asociation Kendall rank correlation test - Kendall’s tau

This is used to estimate a rank-based measure of association. This test may be used if the data do not necessarily come from a bivariate normal distribution.

The correlation coefficient between revenue and spend are 0.7209243 and the p-value < 2.2e-16. Exist correlation

but we know it is not causation

res2

    Kendall's rank correlation tau

data:  df_rev_spend$total_revenue and df_rev_spend$spend
z = 34.731, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.7209243 

We know that they are highly correlated. We can use time series to see the magnitude of the effect of marketing using it as regressor

Visually in the previous plot of revenue, we can see it is stationary, and there is seasonality. We can use auto.arima to calculate the regressor factor

This is a very basic analysis, but it helps to understand in overall how the marketing efforts are creating an effect in the revenue.

1.489071 is the regressor factor

Interpretation of sales increase: For every $1,000 (the unit of ad spend is in thousands USD) increase in advertising spend, the unit sales increase by 1.4891 - It is a positive effect, so ACME probably is using RECAST and the marketing campaign is performaing well

sales_increase

    xreg 
1.489071